home *** CD-ROM | disk | FTP | other *** search
/ Computer Shopper 242 / Issue 242 - April 2008 - DPCS0408DVD.ISO / Software Money Savers / VirtualDub / Source / VirtualDub-1.7.7-src.7z / src / Priss / source / polyphase.cpp < prev   
Encoding:
C/C++ Source or Header  |  2006-03-14  |  32.2 KB  |  1,088 lines
  1. //    Priss (NekoAmp 2.0) - MPEG-1/2 audio decoding library
  2. //    Copyright (C) 2003 Avery Lee
  3. //
  4. //    This program is free software; you can redistribute it and/or modify
  5. //    it under the terms of the GNU General Public License as published by
  6. //    the Free Software Foundation; either version 2 of the License, or
  7. //    (at your option) any later version.
  8. //
  9. //    This program is distributed in the hope that it will be useful,
  10. //    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11. //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12. //    GNU General Public License for more details.
  13. //
  14. //    You should have received a copy of the GNU General Public License
  15. //    along with this program; if not, write to the Free Software
  16. //    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  17.  
  18. #include <string.h>
  19. #include <math.h>
  20. #include <malloc.h>
  21. #include <vd2/system/cpuaccel.h>
  22. #include "polyphase.h"
  23.  
  24. #ifdef _M_AMD64
  25. extern "C" void vdasm_mpegaudio_polyphase_dct4x8(float *f);
  26. extern "C" void vdasm_mpegaudio_polyphase_dctinputbutterflies(float *out, const float *in);
  27. extern "C" void vdasm_mpegaudio_polyphase_matrixout_stereo(const float (*pSrc)[16], const float *pWinFwd, const float *pWinRev, int inc, const uint32 *pSampleInv, const sint16 *pDst, const float (*pSrcFinal)[16], const uint32 *pFinalMask);
  28. #endif
  29.  
  30. extern "C" {
  31.     extern const float __declspec(align(16)) leecoef1[16]={        // 1/(2 cos (pi/64*(2i+1)))
  32.         // first eight (forward)
  33.         0.50060299823520f,
  34.         0.50547095989754f,
  35.         0.51544730992262f,
  36.         0.53104259108978f,
  37.         0.55310389603444f,
  38.         0.58293496820613f,
  39.         0.62250412303566f,
  40.         0.67480834145501f,
  41.  
  42.         // last eight (backward)
  43.         10.19000812354803f,
  44.         3.40760841846872f,
  45.         2.05778100995341f,
  46.         1.48416461631417f,
  47.         1.16943993343288f,
  48.         0.97256823786196f,
  49.         0.83934964541553f,
  50.         0.74453627100230f,
  51.     };
  52.  
  53.     extern const float __declspec(align(16)) leecoef2[8]={        // 1/(2 cos (pi/32*(2i+1)))
  54.         0.50241928618816f,
  55.         0.52249861493969f,
  56.         0.56694403481636f,
  57.         0.64682178335999f,
  58.         0.78815462345125f,
  59.         1.06067768599035f,
  60.         1.72244709823833f,
  61.         5.10114861868916f,
  62.     };
  63. }
  64.  
  65. namespace {
  66.     // Dewindowing coefficients (synthesis filter).
  67.     //
  68.     // There are two significant features of ISO 11172-3 table B.3:
  69.     //
  70.     //   1) All the values are multiples of 1/65536.
  71.     //   2) The FIR filter is antisymmetric about D[256].
  72.     //
  73.     // We store the first half of the coefficients scaled up by 64K.
  74.  
  75.     static const int sFilter[17][16] = {
  76.     {  0, -29, 213, -459, 2037, -5153,  6574,-37489,75038,37489, 6574, 5153, 2037, 459,213,29},
  77.     { -1, -31, 218, -519, 2000, -5517,  5959,-39336,74992,35640, 7134, 4788, 2063, 401,208,26},
  78.     { -1, -35, 222, -581, 1952, -5879,  5288,-41176,74856,33791, 7640, 4425, 2080, 347,202,24},
  79.     { -1, -38, 225, -645, 1893, -6237,  4561,-43006,74630,31947, 8092, 4063, 2087, 294,196,21},
  80.     { -1, -41, 227, -711, 1822, -6589,  3776,-44821,74313,30112, 8492, 3705, 2085, 244,190,19},
  81.     { -1, -45, 228, -779, 1739, -6935,  2935,-46617,73908,28289, 8840, 3351, 2075, 197,183,17},
  82.     { -1, -49, 228, -848, 1644, -7271,  2037,-48390,73415,26482, 9139, 3004, 2057, 153,176,16},
  83.     { -2, -53, 227, -919, 1535, -7597,  1082,-50137,72835,24694, 9389, 2663, 2032, 111,169,14},
  84.     { -2, -58, 224, -991, 1414, -7910,    70,-51853,72169,22929, 9592, 2330, 2001,  72,161,13},
  85.     { -2, -63, 221,-1064, 1280, -8209,  -998,-53534,71420,21189, 9750, 2006, 1962,  36,154,11},
  86.     { -2, -68, 215,-1137, 1131, -8491, -2122,-55178,70590,19478, 9863, 1692, 1919,   2,147,10},
  87.     { -3, -73, 208,-1210,  970, -8755, -3300,-56778,69679,17799, 9935, 1388, 1870, -29,139, 9},
  88.     { -3, -79, 200,-1283,  794, -8998, -4533,-58333,68692,16155, 9966, 1095, 1817, -57,132, 8},
  89.     { -4, -85, 189,-1356,  605, -9219, -5818,-59838,67629,14548, 9959,  814, 1759, -83,125, 7},
  90.     { -4, -91, 177,-1428,  402, -9416, -7154,-61289,66494,12980, 9916,  545, 1698,-106,117, 7},
  91.     { -5, -97, 163,-1498,  185, -9585, -8540,-62684,65290,11455, 9838,  288, 1634,-127,111, 6},
  92.     { -5,-104, 146,-1567,  -45, -9727, -9975,-64019,64019, 9975, 9727,   45, 1567,-146,104, 5},
  93.     };
  94. }
  95.  
  96. void *VDMPEGAudioPolyphaseFilter::operator new(size_t s) {
  97.     return _aligned_malloc(s, 32);
  98. }
  99.  
  100. void VDMPEGAudioPolyphaseFilter::operator delete(void *p) {
  101.     _aligned_free(p);
  102. }
  103.  
  104. VDMPEGAudioPolyphaseFilter *VDMPEGAudioPolyphaseFilter::Create() {
  105.     switch(GetOptMode()) {
  106.     case kOptModeSSE:
  107.         return new VDMPEGAudioPolyphaseFilterSSE;
  108.     default:
  109.         return new VDMPEGAudioPolyphaseFilterFPU;
  110.     }
  111. }
  112.  
  113. VDMPEGAudioPolyphaseFilter::VDMPEGAudioPolyphaseFilter() {
  114.     unsigned i, j;
  115.  
  116.     for(i=0; i<17; ++i) {
  117.         for(j=0; j<16; ++j)
  118.             mFilter[i][j] = mFilter[i][j+16] = sFilter[i][j] * 0.5f; // * (1.0 / 65536.0);
  119.     }
  120.  
  121.     Reset();
  122. }
  123.  
  124. VDMPEGAudioPolyphaseFilter::OptMode VDMPEGAudioPolyphaseFilter::GetOptMode() {
  125.     long optflags = CPUGetEnabledExtensions();
  126.  
  127.     if (optflags & CPUF_SUPPORTS_SSE)
  128.         return kOptModeSSE;
  129.     else
  130.         return kOptModeFPU;
  131. }
  132.  
  133. void VDMPEGAudioPolyphaseFilter::Reset() {
  134.     mWindowPos = 0;
  135.     memset(&mWindow, 0, sizeof mWindow);
  136. }
  137.  
  138. // Matrixing phase
  139. //
  140. // The matrixing phase of 11172-3 is a 32-point DCT transform that produces
  141. // 64 output values.  As there are only 32 input values, 32 of the output
  142. // values are redundant and need not be generated.  To compute the DCT32,
  143. // we apply two iterations of Byeong Gi Lee's decomposition to break it down
  144. // into four 8-point DCTs, then apply the 8-point AAN DCT.
  145. //
  146. // The Lee decomposition consists of five steps:
  147. //
  148. // 1) Butterfly (x+y, x-y) between low and high halves to create even and
  149. //    odd parts.
  150. // 2) Multiply odd part by 1/(2 cos(pi/2N*(2i+1)))
  151. // 3) DCT on both parts.
  152. // 4) Shift-and-add elements in odd half.
  153. // 5) Interleave.
  154. //
  155. // Lee can be used to break down the 8-point DCT as well in 29a12m, but
  156. // it requires three layers of multiplications as well as a high number
  157. // of butterflies, whereas LLM and AAN require less data movement. We do
  158. // not have a quantization stage so AAN is 29a13m instead of 29a5m, but
  159. // LLM has a (2 sqrt 2) scaling factor makes it also 29a13m.  The LLM
  160. // has more rotators which are expensive latency-wise, so we use AAN.
  161.  
  162. void VDMPEGAudioPolyphaseFilter::DCT4x8(float *x) {
  163.     static const float cosvals[8]={        // cos(x*pi/16)
  164.         1.f,
  165.         0.98078528040323043f,
  166.         0.92387953251128674f,
  167.         0.83146961230254524f,
  168.         0.70710678118654757f,
  169.         0.55557023301960229f,
  170.         0.38268343236508984f,
  171.         0.19509032201612833f,
  172.     };
  173.  
  174.     // This is essentially the AAN algorithm, except that leave the rotator
  175.     // as 2a4m, giving us a total of 28a6m.  The operations have been
  176.     // finessed a bit as well to reduce data movement.
  177.  
  178.     // perform B3
  179.  
  180.     float (*const s)[4] = (float (*)[4])x;
  181.     
  182.     for(unsigned i=0; i<4; ++i) {
  183.         // perform B3 (8a)
  184.         const float b3[8] = {
  185.             s[0][i]+s[7][i],
  186.             s[1][i]+s[6][i],
  187.             s[2][i]+s[5][i],
  188.             s[3][i]+s[4][i],
  189.             s[0][i]-s[7][i],
  190.             s[1][i]-s[6][i],
  191.             s[2][i]-s[5][i],
  192.             s[3][i]-s[4][i],
  193.         };
  194.  
  195.         // perform B2, part of ~B1 (7a)
  196.  
  197.         const float b2[8] = {
  198.             b3[0]+b3[3],
  199.             b3[1]+b3[2],
  200.             b3[0]-b3[3],
  201.             b3[1]-b3[2],
  202.             b3[4]-b3[6],
  203.             b3[5]+b3[7],
  204.             b3[6],
  205.             b3[7]-b3[4],
  206.         };
  207.  
  208.         // perform part of ~B1, M (5a6m)
  209.  
  210.         const float m[8] = {
  211.             b2[0]+b2[1],
  212.             b2[0]-b2[1],
  213.             (b2[2]-b2[3])*cosvals[4],
  214.             b2[3],
  215.             b2[4]*cosvals[6] + b2[5]*cosvals[2],
  216.             b2[4]*cosvals[2] - b2[5]*cosvals[6],
  217.             b2[6],
  218.             b2[7]*cosvals[4],
  219.         };
  220.  
  221.         // perform R1 (pass 1) (4a)
  222.  
  223.         const float r1a[8]={
  224.             m[0],
  225.             m[1],
  226.             m[2]+m[3],
  227.             m[2]-m[3],
  228.             m[4],
  229.             m[5],
  230.             m[6]+m[7],
  231.             m[6]-m[7],
  232.         };
  233.  
  234.         // perform R1 (pass 2) (4a)
  235.  
  236.         const float r1[8]={
  237.             r1a[0],
  238.             r1a[1],
  239.             r1a[2],
  240.             r1a[3],
  241.             r1a[7]+r1a[4],
  242.             r1a[6]+r1a[5],
  243.             r1a[6]-r1a[5],
  244.             r1a[7]-r1a[4],
  245.         };
  246.  
  247.         // perform D
  248.  
  249.         static const float d[8]={        // [1 .5 .5 -.5 .5 .5 .5 .5] ./ cos(pi*[0 5 6 1 4 7 2 3]/16)
  250.             1.00000000000000f,
  251.             0.89997622313642f,
  252.             1.30656296487638f,
  253.             -0.50979557910416f,
  254.             0.70710678118655f,
  255.             2.56291544774151f,
  256.             0.54119610014620f,
  257.             0.60134488693505f,
  258.         };
  259.  
  260.         // perform P
  261.  
  262.         s[0][i] = r1[0] * d[0];
  263.         s[1][i] = r1[4] * d[1];
  264.         s[2][i] = r1[2] * d[2];
  265.         s[3][i] = r1[6] * d[3];
  266.         s[4][i] = r1[1] * d[4];
  267.         s[5][i] = r1[5] * d[5];
  268.         s[6][i] = r1[3] * d[6];
  269.         s[7][i] = r1[7] * d[7];
  270.     }
  271. }
  272.  
  273. #ifdef _M_AMD64
  274.     void VDMPEGAudioPolyphaseFilterSSE::DCT4x8(float *x) {
  275.         vdasm_mpegaudio_polyphase_dct4x8(x);
  276.     }
  277. #else
  278.     void __declspec(naked) VDMPEGAudioPolyphaseFilterSSE::DCT4x8(float *x) {
  279.         static const __declspec(align(16)) float c2[4]={0.92387953251128674f,0.92387953251128674f,0.92387953251128674f,0.92387953251128674f};
  280.         static const __declspec(align(16)) float c4[4]={0.70710678118654757f,0.70710678118654757f,0.70710678118654757f,0.70710678118654757f};
  281.         static const __declspec(align(16)) float c6[4]={0.38268343236508984f,0.38268343236508984f,0.38268343236508984f,0.38268343236508984f};
  282.  
  283.         static const __declspec(align(16)) float d[8][4]={        // [1 .5 .5 -.5 .5 .5 .5 .5] ./ cos(pi*[0 5 6 1 4 7 2 3]/16)
  284.     #define TE(x) x,x,x,x
  285.             TE(1.00000000000000f),
  286.             TE(0.89997622313642f),
  287.             TE(1.30656296487638f),
  288.             TE(-0.50979557910416f),
  289.             TE(0.70710678118655f),
  290.             TE(2.56291544774151f),
  291.             TE(0.54119610014620f),
  292.             TE(0.60134488693505f),
  293.     #undef TE
  294.         };
  295.  
  296.         // See the FPU version to get an idea of the flow of this AAN
  297.         // implementation.  Note that we do all four DCTs in parallel!
  298.  
  299.         __asm {
  300.             mov        ecx, [esp+4]
  301.  
  302.             ;even part - B3 (4a)
  303.             movaps    xmm0, [ecx+0*16]        ;xmm0 = s[0]
  304.             movaps    xmm1, [ecx+1*16]        ;xmm1 = s[1]
  305.             movaps    xmm2, [ecx+2*16]        ;xmm2 = s[2]
  306.             movaps    xmm3, [ecx+3*16]        ;xmm3 = s[3]
  307.             addps    xmm0, [ecx+7*16]        ;xmm0 = s[0]+s[7]
  308.             addps    xmm1, [ecx+6*16]        ;xmm1 = s[1]+s[6]
  309.             addps    xmm2, [ecx+5*16]        ;xmm2 = s[2]+s[5]
  310.             addps    xmm3, [ecx+4*16]        ;xmm3 = s[3]+s[4]
  311.  
  312.             ;even part - B2/~B1a (4a)
  313.             movaps    xmm4, xmm0
  314.             addps    xmm0, xmm3                ;xmm0 = b2[0] = b3[0]+b3[3]
  315.             movaps    xmm5, xmm1
  316.             addps    xmm1, xmm2                ;xmm1 = b2[1] = b3[1]+b2[2]
  317.             subps    xmm4, xmm3                ;xmm4 = b2[2] = b3[0]-b3[3]
  318.             subps    xmm5, xmm2                ;xmm5 = b2[3] = b3[1]-b3[2]
  319.  
  320.             ;even part - ~B1b/M (3a1m)
  321.             movaps    xmm2, xmm0
  322.             subps    xmm4, xmm5
  323.             addps    xmm0, xmm1                ;xmm0 = m[0] = b2[0] + b2[1]
  324.             mulps    xmm4, c4                ;xmm4 = m[2] = (b2[2] - b2[3])*c4
  325.             subps    xmm2, xmm1                ;xmm2 = m[1] = b2[0] - b2[1]
  326.  
  327.             ;even part - R1 (2a)
  328.             movaps    xmm3, xmm4
  329.             subps    xmm4, xmm5                ;xmm4 = r1[3] = m[2]-m[3]
  330.             addps    xmm3, xmm5                ;xmm3 = r1[2] = m[2]+m[3]
  331.  
  332.             ;even part - d (4m)
  333.             mulps    xmm0, [d+0*16]            ;xmm0 = out[0] = r1[0]*d[0]
  334.             mulps    xmm2, [d+4*16]            ;xmm2 = out[4] = r1[1]*d[4]
  335.             mulps    xmm3, [d+2*16]            ;xmm3 = out[2] = r1[2]*d[2]
  336.             mulps    xmm4, [d+6*16]            ;xmm4 = out[6] = r1[3]*d[6]
  337.  
  338.             ;odd part - B3 (4a)
  339.             movaps    xmm1, [ecx+0*16]
  340.             movaps    xmm5, [ecx+1*16]
  341.             movaps    xmm6, [ecx+2*16]
  342.             movaps    xmm7, [ecx+3*16]
  343.             subps    xmm1, [ecx+7*16]        ;xmm1 = b3[4] = s[0]-s[7]
  344.             subps    xmm5, [ecx+6*16]        ;xmm5 = b3[5] = s[1]-s[6]
  345.             subps    xmm6, [ecx+5*16]        ;xmm6 = b3[6] = s[2]-s[5]
  346.             subps    xmm7, [ecx+4*16]        ;xmm7 = b3[7] = s[3]-s[4]
  347.  
  348.             ;even part - writeout
  349.             movaps    [ecx+0*16], xmm0
  350.             movaps    [ecx+4*16], xmm2
  351.             movaps    [ecx+2*16], xmm3
  352.             movaps    [ecx+6*16], xmm4
  353.  
  354.             ;odd part - B2/~B1a (3a)
  355.             addps    xmm5, xmm7                ;xmm5 = b2[5] = b3[5]+b3[7]
  356.             subps    xmm7, xmm1                ;xmm7 = b2[7] = b3[7]-b3[4]
  357.             subps    xmm1, xmm6                ;xmm1 = b2[4] = b3[4]-b3[6]
  358.  
  359.             ;odd part - ~B1b/M (2a5m)
  360.             movaps    xmm0, xmm1
  361.             mulps    xmm7, c4                ;xmm7 = m[7] = c4*b2[7]
  362.             movaps    xmm2, xmm5
  363.             mulps    xmm0, c6
  364.             mulps    xmm1, c2
  365.             mulps    xmm2, c2
  366.             mulps    xmm5, c6
  367.             addps    xmm0, xmm2                ;xmm0 = m[4] = c6*b2[4] + c2*b2[5]
  368.             subps    xmm1, xmm5                ;xmm1 = m[5] = c2*b2[4] - c6*b2[5]
  369.  
  370.             ;odd part - R1a (2a)
  371.             movaps    xmm5, xmm6
  372.             addps    xmm6, xmm7                ;xmm6 = r1a[6] = m[6]+m[7]
  373.             subps    xmm5, xmm7                ;xmm5 = r1a[7] = m[6]-m[7]
  374.  
  375.             ;odd part - R1b (4a)
  376.             movaps    xmm3, xmm5
  377.             movaps    xmm4, xmm6
  378.             subps    xmm5, xmm0                ;xmm5 = r1b[7] = r1a[7]-r1a[4]
  379.             subps    xmm6, xmm1                ;xmm6 = r1b[6] = r1a[6]-r1a[5]
  380.             addps    xmm4, xmm1                ;xmm4 = r1b[5] = r1a[6]+r1a[5]
  381.             addps    xmm3, xmm0                ;xmm3 = r1b[4] = r1a[7]+r1a[4]
  382.  
  383.             ;odd part - D (4a)
  384.             mulps    xmm3, [d+1*16]
  385.             mulps    xmm4, [d+5*16]
  386.             mulps    xmm6, [d+3*16]
  387.             mulps    xmm5, [d+7*16]
  388.  
  389.             ;odd part - writeout
  390.             movaps    [ecx+1*16], xmm3
  391.             movaps    [ecx+5*16], xmm4
  392.             movaps    [ecx+3*16], xmm6
  393.             movaps    [ecx+7*16], xmm5
  394.  
  395.             ret        4
  396.         }
  397.     }
  398. #endif
  399.  
  400. namespace {
  401. #ifdef _M_AMD64
  402.     static void __stdcall DCTInputButterfliesSSE(float x[32], const float in[32]) {
  403.         vdasm_mpegaudio_polyphase_dctinputbutterflies(x, in);
  404.     }
  405. #else
  406.     static void __declspec(naked) __stdcall DCTInputButterfliesSSE(float x[32], const float in[32]) {
  407.         __asm {
  408.             push    ebx
  409.             mov        eax, [esp+8+4]
  410.             mov        ebx, [esp+4+4]
  411.             xor        ecx, ecx
  412.             mov        edx, 48
  413. xloop:
  414.             movups    xmm0, [eax+ecx]            ;xmm0 = in[i]
  415.             movups    xmm1, [eax+edx]            ;xmm1 = in[15-i]
  416.             movups    xmm2, [eax+ecx+64]        ;xmm2 = in[i+16]
  417.             movups    xmm3, [eax+edx+64]        ;xmm3 = in[31-i]
  418.             shufps    xmm1, xmm1, 00011011b
  419.             shufps    xmm3, xmm3, 00011011b
  420.  
  421.             ;butterfly for first decomposition
  422.             movaps    xmm4, xmm0
  423.             movaps    xmm5, xmm1
  424.             addps    xmm0, xmm3                ;xmm0 = y0 = x0+x3
  425.             addps    xmm1, xmm2                ;xmm1 = y1 = x1+x2
  426.             subps    xmm4, xmm3                ;xmm4 = y2 = x0-x3
  427.             subps    xmm5, xmm2                ;xmm5 = y3 = x1-x2
  428.             mulps    xmm4, [leecoef1+ecx]
  429.             mulps    xmm5, [leecoef1+ecx+32]
  430.  
  431.             ;butterfly for second decomposition
  432.             movaps    xmm2, xmm0
  433.             movaps    xmm3, xmm4
  434.             addps    xmm0, xmm1                ;xmm0 = z0 = y0+y1
  435.             subps    xmm2, xmm1                ;xmm2 = z1 = y0-y1
  436.             addps    xmm3, xmm5                ;xmm3 = z2 = y2+y3
  437.             subps    xmm4, xmm5                ;xmm4 = z3 = y2-y3
  438.             mulps    xmm2, [leecoef2+ecx]
  439.             mulps    xmm4, [leecoef2+ecx]
  440.  
  441.             ;interleave in 0-2-1-3 order
  442.             movaps        xmm1, xmm0
  443.             unpcklps    xmm0, xmm3            ;xmm0 = z2B | z0B | z2A | z0A
  444.             unpckhps    xmm1, xmm3            ;xmm1 = z2D | z0D | z2C | z0C
  445.             movaps        xmm6, xmm2
  446.             unpcklps    xmm2, xmm4            ;xmm2 = z3B | z1B | z3A | z1A
  447.             unpckhps    xmm6, xmm4            ;xmm6 = z3D | z1D | z3C | z1C
  448.  
  449.             movlps    [ebx   ], xmm0
  450.             movlps    [ebx+ 8], xmm2
  451.             movhps    [ebx+16], xmm0
  452.             movhps    [ebx+24], xmm2
  453.             movlps    [ebx+32], xmm1
  454.             movlps    [ebx+40], xmm6
  455.             movhps    [ebx+48], xmm1
  456.             movhps    [ebx+56], xmm6
  457.  
  458.             add        ebx, 64
  459.             add        ecx, 16
  460.             sub        edx, 16
  461.             cmp        ecx, edx
  462.             jb        xloop
  463.  
  464.             pop        ebx
  465.             ret        8
  466.         }
  467.     }
  468. #endif
  469. }
  470.  
  471. void VDMPEGAudioPolyphaseFilterSSE::DCTInputButterflies(float x[32], const float in[32]) {
  472.     DCTInputButterfliesSSE(x, in);
  473. }
  474.  
  475. void VDMPEGAudioPolyphaseFilter::DCTInputButterflies(float x[32], const float in[32]) {
  476.     for(unsigned i=0; i<8; ++i) {
  477.         const float x0 = in[i];
  478.         const float x1 = in[15-i];
  479.         const float x2 = in[i+16];
  480.         const float x3 = in[31-i];
  481.  
  482.         // butterfly for first decomposition
  483.         const float y0 =  x0+x3;
  484.         const float y1 =  x1+x2;
  485.         const float y2 = (x0-x3)*leecoef1[i];
  486.         const float y3 = (x1-x2)*leecoef1[i+8];        // last 8 are reflected in table
  487.  
  488.         // butterfly for second decomposition
  489.         const float z0 =  y0+y1;
  490.         const float z1 = (y0-y1)*leecoef2[i];
  491.         const float z2 =  y2+y3;
  492.         const float z3 = (y2-y3)*leecoef2[i];
  493.  
  494.         // store with output reordering
  495.         x[i*4+0] = z0;
  496.         x[i*4+1] = z2;
  497.         x[i*4+2] = z1;
  498.         x[i*4+3] = z3;
  499.     }
  500. }
  501.  
  502. void VDMPEGAudioPolyphaseFilter::Matrix(const float *in, bool stereo, int ch) {
  503.     unsigned i;
  504.     float __declspec(align(16)) x[32];
  505.  
  506.     // do input butterflies to split down to four 8-point DCTs
  507.     DCTInputButterflies(x, in);
  508.  
  509.     // do four 8-point AAN DCTs
  510.     DCT4x8(x);
  511.  
  512.     // offset addition for odd 8-point matrices
  513.     x[ 2] += x[ 6];
  514.     x[ 3] += x[ 7];
  515.  
  516.     x[ 6] += x[10];
  517.     x[ 7] += x[11];
  518.  
  519.     x[10] += x[14];
  520.     x[11] += x[15];
  521.  
  522.     x[14] += x[18];
  523.     x[15] += x[19];
  524.  
  525.     x[18] += x[22];
  526.     x[19] += x[23];
  527.  
  528.     x[22] += x[26];
  529.     x[23] += x[27];
  530.  
  531.     x[26] += x[30];
  532.     x[27] += x[31];
  533.  
  534.     // offset addition for odd 16-point matrix, as well as reflection of odd
  535.     // sample sets to aid the synthesis filter bank
  536.  
  537.     if (stereo) {
  538.         float (*out)[2][16] = (float (*)[2][16])&mWindow.stereo[0][ch][mWindowPos];
  539.  
  540.         if (mWindowPos & 1) {
  541.             out[0 ][0][0] = x[0];
  542.             out[31][0][0] = x[1] + x[3];
  543.  
  544.             for (i = 2; i < 30; i+=2) {
  545.                 out[32-i][0][0] = x[i];
  546.                 out[31-i][0][0] = x[i+1] + x[i+3];
  547.             }
  548.  
  549.             out[2][0][0] = x[30];
  550.             out[1][0][0] = x[31];
  551.         } else {
  552.             for (i = 0; i < 30; i+=2) {
  553.                 out[i+0][0][0] = x[i];
  554.                 out[i+1][0][0] = x[i+1] + x[i+3];
  555.             }
  556.  
  557.             out[30][0][0] = x[30];
  558.             out[31][0][0] = x[31];
  559.         }
  560.     } else {
  561.         float (*out)[16] = (float (*)[16])&mWindow.mono[0][mWindowPos];
  562.  
  563.         if (mWindowPos & 1) {
  564.             out[0 ][0] = x[0];
  565.             out[31][0] = x[1] + x[3];
  566.  
  567.             for (i = 2; i < 30; i+=2) {
  568.                 out[32-i][0] = x[i];
  569.                 out[31-i][0] = x[i+1] + x[i+3];
  570.             }
  571.  
  572.             out[2][0] = x[30];
  573.             out[1][0] = x[31];
  574.         } else {
  575.             for (i = 0; i < 30; i+=2) {
  576.                 out[i+0][0] = x[i];
  577.                 out[i+1][0] = x[i+1] + x[i+3];
  578.             }
  579.  
  580.             out[30][0] = x[30];
  581.             out[31][0] = x[31];
  582.         }
  583.     }
  584. }
  585.  
  586. // Subband synthesis filter strategy
  587. //
  588. // Since we only stored half of the 64-point DCT output in the synthesis
  589. // window, we must emulate the 32>64 conversion. From Intel's AP-533:
  590. //
  591. //        x[ 0..15] =  x'[16..31]
  592. //      x[16]     = 0
  593. //      x[17..31] = -x'[31..17]
  594. //      x[32..48] = -x'[16.. 0]
  595. //        x[48..63] = -x'[ 0..15]
  596. //
  597. // Even offsets use x[0..31], while odd offsets use x[32..63].  There is
  598. // symmetry in the sample generation:
  599. //
  600. // even x[ 0]     =  x'[16]
  601. //        x[ 1..15] =  x'[16 + 1..15]
  602. //      x[16]     = 0
  603. //      x[17..31] = -x'[16 + 15..1]
  604. //
  605. // odd  x[32]     = -x'[16]
  606. //      x[33..47] = -x'[16 - (1..15)]
  607. //      x[48]     = -x'[ 0]
  608. //        x[49..63] = -x'[16 - (15..1)]
  609. //
  610. // So we can split the reconstruction into three phases: s[0], s[16], and
  611. // a phase that generates s[1..15] and s[31..17] at the same time, in
  612. // opposite directions.  As it turns out, the subband synthesis window is
  613. // antisymmetric around D[256], meaning that row 17 is the negative
  614. // reverse of row 15, row 18 is the negative reverse of row 14, etc. and
  615. // we can reuse window rows in reverse as well.  We now have only 17
  616. // window rows to deal with.
  617. //
  618. // Getting rid of the sliding window is easy: use a circular array and
  619. // double up all the window rows so we can choose a contiguous array of
  620. // 16 window elements for any offset of the sliding window.  Next.
  621. //
  622. // That leaves the question of how to deal with the annoying flip/flop
  623. // of the reconstruction filter between the first and second halves of
  624. // the matrixing stage output.  The V[] matrix looks like this in
  625. // the standard:
  626. //
  627. //        A1
  628. //            A2
  629. //        B1
  630. //            B2
  631. //        C1
  632. //            C2
  633. //        D1
  634. //            D2
  635. //
  636. // On the first run we would process A1-B2-C1-D2, on the second run
  637. // B1-C2-D1-E2, and so on.  However, if we swap the halves on every
  638. // other block:
  639. //
  640. //        A1
  641. //            A2
  642. //        B2
  643. //            B1
  644. //        C1
  645. //            C2
  646. //        D2
  647. //            D1
  648. //
  649. // Then the hopping is easy -- we just process all primary halves on even
  650. // passes and all secondary halves on odd passes.
  651. //
  652. // Finally, to eliminate the redundant 32 elements per set: notice that
  653. // our traversal over odd elements is mirrored from even ones -- the first
  654. // 16 samples draw from x[16..31] on even sets and x[16..1] on odd sets,
  655. // mirrored around x[16].  So we just mirror every other set on write and
  656. // mirror direction for every other set on read.  We will unfortunately
  657. // have to have two loops for sliding window offset polarity since
  658. // x[0..15] and x[33..48] are opposite sign but x[17..31] and x[48..63]
  659. // are not.
  660.  
  661. namespace {
  662. #if 1
  663.     int clip(float f) {
  664.         static const float bias = 8388608.0f + 4194304.0f + 32768.0f; // 2^23 + 2^22
  665.         float tmp = f + bias;
  666.         int v = reinterpret_cast<sint32&>(tmp) - 0x4b400000;
  667.  
  668.         if ((unsigned)v >= 65536)
  669.             v = (~v >> 31)&0xffff;
  670.  
  671.         return v - 0x8000;
  672.     }
  673. #else
  674.     int clip(float f) {
  675.         int v = (int)f;
  676.  
  677.         if (v < -32768)
  678.             v = -32768;
  679.         else if (v > 32767)
  680.             v = 32767;
  681.  
  682.         return v;
  683.     }
  684. #endif
  685. }
  686.  
  687. void VDMPEGAudioPolyphaseFilter::SynthesizeMono(sint16 *dst) {
  688.     bool odd = (mWindowPos & 1) != 0;
  689.     int sample;
  690.  
  691.     // first sample is special case
  692.     {
  693.         const float *window_fwd = &mFilter[0][(-(int)mWindowPos)&15];
  694.         const float *src = &mWindow.mono[16][0];
  695.         float left = 0;
  696.  
  697.         for(int i=0; i<16; i+=2) {
  698.             left  += src[i+0] * window_fwd[i];
  699.             left  -= src[i+1] * window_fwd[i+1];
  700.         }
  701.  
  702.         if (odd) {
  703.             left = -left;
  704.         }
  705.  
  706.         dst[ 0] = clip(left);
  707.     }
  708.  
  709.     // reflected samples
  710.     for(sample=1; sample<16; ++sample) {
  711.         const float *window_fwd = &mFilter[sample][(-(int)mWindowPos)&15];
  712.         const float *window_rev = &mFilter[sample][mWindowPos&15];
  713.         const float *src = mWindow.mono[odd ? 16-sample : 16+sample];
  714.         float left1 = 0, left2 = 0;
  715.  
  716.         for(int i=0; i<16; i+=2) {
  717.             left1  += src[i+0] * window_fwd[i];
  718.             left1  -= src[i+1] * window_fwd[i+1];
  719.  
  720.             left2  += src[i+0] * window_rev[15-i];
  721.             left2  += src[i+1] * window_rev[14-i];
  722.         }
  723.  
  724.         if (odd) {
  725.             left1 = -left1;
  726.         }
  727.  
  728.         dst[ 0+sample] = clip(left1);
  729.         dst[32-sample] = clip(left2);
  730.     }
  731.  
  732.     // center sample is special case
  733.     {
  734.         const float *window_fwd = &mFilter[16][(-(int)mWindowPos)&15];
  735.         const float *src = mWindow.mono[0];
  736.         float left = 0;
  737.  
  738.         for(int i=!odd; i<16; i+=2) {
  739.             left  -= src[i] * window_fwd[i];
  740.         }
  741.  
  742.         dst[16] = clip(left);
  743.     }
  744. }
  745.  
  746. void VDMPEGAudioPolyphaseFilter::SynthesizeStereo(sint16 *dst) {
  747.     bool odd = (mWindowPos & 1) != 0;
  748.     int sample;
  749.  
  750.     // first sample is special case
  751.     {
  752.         const float *window_fwd = &mFilter[0][(-(int)mWindowPos)&15];
  753.         const float (*src)[16] = &mWindow.stereo[16][0];
  754.         float left = 0, right = 0;
  755.  
  756.         for(int i=0; i<16; i+=2) {
  757.             left  += src[0][i+0] * window_fwd[i];
  758.             right += src[1][i+0] * window_fwd[i];
  759.             left  -= src[0][i+1] * window_fwd[i+1];
  760.             right -= src[1][i+1] * window_fwd[i+1];
  761.         }
  762.  
  763.         if (odd) {
  764.             left = -left;
  765.             right = -right;
  766.         }
  767.  
  768.         dst[ 0] = clip(left);
  769.         dst[ 1] = clip(right);
  770.     }
  771.  
  772.     // reflected samples
  773.     for(sample=1; sample<16; ++sample) {
  774.         const float *window_fwd = &mFilter[sample][(-(int)mWindowPos)&15];
  775.         const float *window_rev = &mFilter[sample][mWindowPos&15];
  776.         const float (*src)[16] = mWindow.stereo[odd ? 16-sample : 16+sample];
  777.         float left1 = 0, right1 = 0, left2 = 0, right2 = 0;
  778.  
  779.         for(int i=0; i<16; i+=2) {
  780.             left1  += src[0][i+0] * window_fwd[i];
  781.             right1 += src[1][i+0] * window_fwd[i];
  782.             left1  -= src[0][i+1] * window_fwd[i+1];
  783.             right1 -= src[1][i+1] * window_fwd[i+1];
  784.  
  785.             left2  += src[0][i+0] * window_rev[15-i];
  786.             right2 += src[1][i+0] * window_rev[15-i];
  787.             left2  += src[0][i+1] * window_rev[14-i];
  788.             right2 += src[1][i+1] * window_rev[14-i];
  789.         }
  790.  
  791.         if (odd) {
  792.             left1 = -left1;
  793.             right1 = -right1;
  794.         }
  795.  
  796.         dst[ 0+sample*2] = clip(left1);
  797.         dst[ 1+sample*2] = clip(right1);
  798.         dst[64-sample*2] = clip(left2);
  799.         dst[65-sample*2] = clip(right2);
  800.     }
  801.  
  802.     // center sample is special case
  803.     {
  804.         const float *window_fwd = &mFilter[16][(-(int)mWindowPos)&15];
  805.         const float (*src)[16] = &mWindow.stereo[0][0];
  806.         float left = 0, right = 0;
  807.  
  808.         for(int i=!odd; i<16; i+=2) {
  809.             left  -= src[0][i] * window_fwd[i];
  810.             right -= src[1][i] * window_fwd[i];
  811.         }
  812.  
  813.         dst[32] = clip(left);
  814.         dst[33] = clip(right);
  815.     }
  816. }
  817.  
  818. namespace {
  819. #ifndef _M_AMD64
  820.     void __declspec(naked) __stdcall ComputeSamplesStereoSSE(const float (*pSrc)[16], const float *pWinFwd, const float *pWinRev, int inc, const uint32 *pSampleInv, const sint16 *pDst, const float (*pSrcFinal)[16], const uint32 *pFinalMask) {
  821.         static const __declspec(align(16)) uint32 invother[4]={ 0, 0x80000000, 0, 0x80000000 };
  822.         static const __declspec(align(16)) uint32 invall[4]={ 0x80000000, 0x80000000, 0x80000000, 0x80000000 };
  823.         __asm {
  824.             push    ebp
  825.             push    edi
  826.             push    esi
  827.             push    ebx
  828.  
  829.             mov        eax, [esp+4+16]            ;eax = pointer to subband samples
  830.             mov        ebx, [esp+20+16]        ;ebx = pointer to sample inversion value
  831.             mov        ecx, [esp+8+16]            ;ecx = pointer to forward window
  832.             mov        edx, [esp+12+16]        ;edx = pointer to reverse window
  833.             mov        esi, [esp+24+16]        ;esi = pointer to first two forward destination samples
  834.             mov        ebp, [esp+16+16]        ;ebp = source increment
  835.             lea        edi, [esi+120]            ;edi = pointer to first two reverse destination samples
  836.  
  837.             ;compute first sample (0)
  838.  
  839.             movaps    xmm5, invother
  840.             movups    xmm0, [ecx]                ;load window samples 0-3
  841.             xorps    xmm0, xmm5                ;toggle signs on odd window samples
  842.             movaps    xmm1, xmm0
  843.             mulps    xmm0, [eax]                ;multiply by left subband samples
  844.             mulps    xmm1, [eax+64]            ;multiply by right subband samples
  845.             movups    xmm2, [ecx+16]            ;load window samples 4-7
  846.             xorps    xmm2, xmm5                ;toggle signs on odd window samples
  847.             movaps    xmm3, xmm2
  848.             mulps    xmm2, [eax+16]            ;multiply by left subband samples
  849.             mulps    xmm3, [eax+80]            ;multiply by right subband samples
  850.             addps    xmm0, xmm2
  851.             addps    xmm1, xmm3
  852.             movups    xmm2, [ecx+32]            ;load window samples 8-11
  853.             xorps    xmm2, xmm5                ;toggle signs on odd window samples
  854.             movaps    xmm3, xmm2
  855.             mulps    xmm2, [eax+32]            ;multiply by left subband samples
  856.             mulps    xmm3, [eax+96]            ;multiply by right subband samples
  857.             addps    xmm0, xmm2
  858.             addps    xmm1, xmm3
  859.             movups    xmm2, [ecx+48]            ;load window samples 12-15
  860.             xorps    xmm2, xmm5                ;toggle signs on odd window samples
  861.             movaps    xmm3, xmm2
  862.             mulps    xmm2, [eax+48]            ;multiply by left subband samples
  863.             mulps    xmm3, [eax+112]            ;multiply by right subband samples
  864.             addps    xmm0, xmm2
  865.             addps    xmm1, xmm3
  866.  
  867.             movaps    xmm2, xmm0                ;xmm2 = l3 | l2 | l1 | l0
  868.             movlhps    xmm0, xmm1                ;xmm0 = r1 | r0 | l1 | l0
  869.             movhlps    xmm1, xmm2                ;xmm1 = r3 | r2 | l3 | l2
  870.             addps    xmm0, xmm1                ;xmm0 = r1+r3 | r0+r2 | l1+l3 | l0+l2
  871.             shufps    xmm0, xmm0, 11011000b    ;xmm0 = r1+r3 | l1+l3 | r0+r2 | l0+l2
  872.             movhlps    xmm3, xmm0                ;xmm3 =   ?   |   ?   | r1+r3 | l1+l3
  873.             movaps    xmm4, [ebx]
  874.             movhlps    xmm4, xmm4
  875.             addps    xmm0, xmm3                ;xmm0 = ? | ? | r | l
  876.             xorps    xmm0, xmm4
  877.             cvtps2pi    mm0, xmm0
  878.             packssdw    mm0, mm0
  879.             movd    [esi-4], mm0
  880.  
  881.             add        ecx, 128
  882.             add        edx, 128
  883.             add        eax, ebp
  884.  
  885.             ;compute reflected samples (1-15, 17-31)
  886. xloop:
  887.             movups    xmm2, [edx+48]
  888.             shufps    xmm2, xmm2, 00011011b    ;xmm2 = reverse window
  889.             movups    xmm3, [ecx]                ;xmm3 = forward window
  890.             xorps    xmm3, invother            ;negate every other sample in forward window
  891.             movaps    xmm0, [eax]                ;xmm0 = left source
  892.             movaps    xmm1, xmm0
  893.             mulps    xmm0, xmm2
  894.             mulps    xmm1, xmm3
  895.             movaps    xmm4, xmm0                ;xmm4 = left forward
  896.             movaps    xmm5, xmm1                ;xmm5 = left reverse
  897.             movaps    xmm0, [eax+64]            ;xmm0 = left source
  898.             movaps    xmm1, xmm0
  899.             mulps    xmm0, xmm2
  900.             mulps    xmm1, xmm3
  901.             movaps    xmm6, xmm0                ;xmm6 = right forward
  902.             movaps    xmm7, xmm1                ;xmm7 = right reverse
  903.  
  904.             movups    xmm2, [edx+32]
  905.             shufps    xmm2, xmm2, 00011011b    ;xmm2 = reverse window
  906.             movups    xmm3, [ecx+16]            ;xmm3 = forward window
  907.             xorps    xmm3, invother            ;negate every other sample in forward window
  908.             movaps    xmm0, [eax+16]            ;xmm0 = left source
  909.             movaps    xmm1, xmm0
  910.             mulps    xmm0, xmm2
  911.             mulps    xmm1, xmm3
  912.             addps    xmm4, xmm0                ;xmm4 += left forward
  913.             addps    xmm5, xmm1                ;xmm5 += left reverse
  914.             movaps    xmm0, [eax+80]            ;xmm0 = left source
  915.             movaps    xmm1, xmm0
  916.             mulps    xmm0, xmm2
  917.             mulps    xmm1, xmm3
  918.             addps    xmm6, xmm0                ;xmm6 += right forward
  919.             addps    xmm7, xmm1                ;xmm7 += right reverse
  920.  
  921.             movups    xmm2, [edx+16]
  922.             shufps    xmm2, xmm2, 00011011b    ;xmm2 = reverse window
  923.             movups    xmm3, [ecx+32]            ;xmm3 = forward window
  924.             xorps    xmm3, invother            ;negate every other sample in forward window
  925.             movaps    xmm0, [eax+32]            ;xmm0 = left source
  926.             movaps    xmm1, xmm0
  927.             mulps    xmm0, xmm2
  928.             mulps    xmm1, xmm3
  929.             addps    xmm4, xmm0                ;xmm4 += left forward
  930.             addps    xmm5, xmm1                ;xmm5 += left reverse
  931.             movaps    xmm0, [eax+96]            ;xmm0 = left source
  932.             movaps    xmm1, xmm0
  933.             mulps    xmm0, xmm2
  934.             mulps    xmm1, xmm3
  935.             addps    xmm6, xmm0                ;xmm6 += right forward
  936.             addps    xmm7, xmm1                ;xmm7 += right reverse
  937.  
  938.             movups    xmm2, [edx]
  939.             shufps    xmm2, xmm2, 00011011b    ;xmm2 = reverse window
  940.             movups    xmm3, [ecx+48]            ;xmm3 = forward window
  941.             xorps    xmm3, invother            ;negate every other sample in forward window
  942.             movaps    xmm0, [eax+48]            ;xmm0 = left source
  943.             movaps    xmm1, xmm0
  944.             mulps    xmm0, xmm2
  945.             mulps    xmm1, xmm3
  946.             addps    xmm4, xmm0                ;xmm4 += left forward
  947.             addps    xmm5, xmm1                ;xmm5 += left reverse
  948.             movaps    xmm0, [eax+112]            ;xmm0 = left source
  949.             movaps    xmm1, xmm0
  950.             mulps    xmm0, xmm2
  951.             mulps    xmm1, xmm3
  952.             addps    xmm6, xmm0                ;xmm6 += right forward
  953.             addps    xmm7, xmm1                ;xmm7 += right reverse
  954.  
  955.             movaps    xmm0, xmm4                ;xmm0 = lf3 | lf2 | lf1 | lf0
  956.             movaps    xmm1, xmm5
  957.             movlhps    xmm0, xmm6                ;xmm0 = rf0 | rf1 | lf1 | lf0
  958.             movlhps    xmm1, xmm7
  959.             movhlps    xmm6, xmm4                ;xmm6 = rf3 | rf2 | lf3 | lf2
  960.             movhlps    xmm7, xmm5
  961.             addps    xmm0, xmm6                ;xmm0 = rf0+rf3 | rf1+rf2 | lf1+lf3 | lf0+lf2
  962.             addps    xmm1, xmm7
  963.             movaps    xmm2, xmm0
  964.             movaps    xmm3, xmm1
  965.             shufps    xmm0, xmm0, 10110001b    ;xmm0 = rf1+rf2 | rf0+rf3 | lf0+lf2 | lf1+lf3
  966.             shufps    xmm1, xmm1, 10110001b
  967.             addps    xmm0, xmm2                ;xmm0 = rf | rf | lf | lf
  968.             addps    xmm1, xmm3                ;xmm1 = rb | rb | lb | lb
  969.             shufps    xmm0, xmm1, 10001000b    ;xmm0 = rf | lf | rb | lb
  970.             xorps    xmm0, [ebx]
  971.             cvtps2pi    mm0, xmm0
  972.             movhlps    xmm0, xmm0
  973.             packssdw    mm0, mm0
  974.             cvtps2pi    mm1, xmm0
  975.             packssdw    mm1, mm1
  976.             movd    [esi], mm1
  977.             movd    [edi], mm0
  978.  
  979.             add        esi,4
  980.             sub        edi,4
  981.             add        eax,ebp
  982.             add        ecx,128
  983.             add        edx,128
  984.             cmp        esi,edi
  985.             jne        xloop
  986.  
  987.             ;do last sample (16)
  988.             mov        eax, [esp+28+16]
  989.             mov        edi, [esp+32+16]
  990.             movaps    xmm5, [edi]                ;load final mask (masks out every other sample)
  991.  
  992.             movups    xmm0, [ecx]                ;load window samples 0-3
  993.             andps    xmm0, xmm5                ;mask out every other sample
  994.             movaps    xmm1, xmm0
  995.             mulps    xmm0, [eax]                ;multiply by left subband samples
  996.             mulps    xmm1, [eax+64]            ;multiply by right subband samples
  997.             movups    xmm2, [ecx+16]            ;load window samples 4-7
  998.             andps    xmm2, xmm5                ;mask out every other sample
  999.             movaps    xmm3, xmm2
  1000.             mulps    xmm2, [eax+16]            ;multiply by left subband samples
  1001.             mulps    xmm3, [eax+80]            ;multiply by right subband samples
  1002.             addps    xmm0, xmm2
  1003.             addps    xmm1, xmm3
  1004.             movups    xmm2, [ecx+32]            ;load window samples 8-11
  1005.             andps    xmm2, xmm5                ;mask out every other sample
  1006.             movaps    xmm3, xmm2
  1007.             mulps    xmm2, [eax+32]            ;multiply by left subband samples
  1008.             mulps    xmm3, [eax+96]            ;multiply by right subband samples
  1009.             addps    xmm0, xmm2
  1010.             addps    xmm1, xmm3
  1011.             movups    xmm2, [ecx+48]            ;load window samples 12-15
  1012.             andps    xmm2, xmm5                ;mask out every other sample
  1013.             movaps    xmm3, xmm2
  1014.             mulps    xmm2, [eax+48]            ;multiply by left subband samples
  1015.             mulps    xmm3, [eax+112]            ;multiply by right subband samples
  1016.             addps    xmm0, xmm2
  1017.             addps    xmm1, xmm3
  1018.  
  1019.             movaps    xmm2, xmm0                ;xmm2 = l3 | l2 | l1 | l0
  1020.             movlhps    xmm0, xmm1                ;xmm0 = r1 | r0 | l1 | l0
  1021.             movhlps    xmm1, xmm2                ;xmm1 = r3 | r2 | l3 | l2
  1022.             addps    xmm0, xmm1                ;xmm0 = r1+r3 | r0+r2 | l1+l3 | l0+l2
  1023.             shufps    xmm0, xmm0, 11011000b    ;xmm0 = r1+r3 | l1+l3 | r0+r2 | l0+l2
  1024.             movhlps    xmm3, xmm0                ;xmm3 =   ?   |   ?   | r1+r3 | l1+l3
  1025.             addps    xmm0, xmm3                ;xmm0 = ? | ? | r | l
  1026.             xorps    xmm0, invall
  1027.             cvtps2pi    mm0, xmm0
  1028.             packssdw    mm0, mm0
  1029.             movd    [esi], mm0
  1030.  
  1031.             emms
  1032.             pop        ebx
  1033.             pop        esi
  1034.             pop        edi
  1035.             pop        ebp
  1036.             ret        32
  1037.         }
  1038.     }
  1039. #endif
  1040. };
  1041.  
  1042. void VDMPEGAudioPolyphaseFilterSSE::SynthesizeStereo(sint16 *dst) {
  1043.     bool odd = (mWindowPos & 1) != 0;
  1044.  
  1045.     // reflected samples
  1046.     static const __declspec(align(16)) uint32 invert[4]={0,0,0x80000000,0x80000000};
  1047.     static const __declspec(align(16)) uint32 noinvert[4]={0,0,0,0};
  1048.     static const __declspec(align(16)) uint32 evenmask[4]={~0,0,~0,0};
  1049.     static const __declspec(align(16)) uint32 oddmask[4]={0,~0,0,~0};
  1050.  
  1051. #ifdef _M_AMD64
  1052.     vdasm_mpegaudio_polyphase_matrixout_stereo(
  1053.                     mWindow.stereo[16],
  1054.                     &mFilter[0][(-mWindowPos)&15],
  1055.                     &mFilter[0][mWindowPos&15],
  1056.                     odd ? -sizeof(mWindow.stereo[0]) : +sizeof(mWindow.stereo[0]),
  1057.                     odd ? invert : noinvert,
  1058.                     &dst[2],
  1059.                     mWindow.stereo[0],
  1060.                     odd ? evenmask : oddmask);
  1061. #else
  1062.     ComputeSamplesStereoSSE(
  1063.                     mWindow.stereo[16],
  1064.                     &mFilter[0][(-(int)mWindowPos)&15],
  1065.                     &mFilter[0][mWindowPos&15],
  1066.                     odd ? -(int)sizeof(mWindow.stereo[0]) : +sizeof(mWindow.stereo[0]),
  1067.                     odd ? invert : noinvert,
  1068.                     &dst[2],
  1069.                     mWindow.stereo[0],
  1070.                     odd ? evenmask : oddmask);
  1071. #endif
  1072. }
  1073.  
  1074. void VDMPEGAudioPolyphaseFilter::Generate(const float left[32], const float right[32], sint16 *dst) {
  1075.     if (right) {
  1076.         Matrix(left, true, 0);
  1077.         Matrix(right, true, 1);
  1078.  
  1079.         SynthesizeStereo(dst);
  1080.     } else {
  1081.         Matrix(left, false, 0);
  1082.  
  1083.         SynthesizeMono(dst);
  1084.     }
  1085.  
  1086.     mWindowPos = (mWindowPos-1) & 15;
  1087. }
  1088.